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Abstract 

We discuss the possibility of solving the inverse problem of the cosmic-ray elec- 
tron/positron from decaying/annihilating dark matter, and show simple analytic 
formulae to reconstruct the source spectrum of the electron/positron from the ob- 
served flux. We also illustrate our approach by applying the obtained formula to 
the just released Fermi data as well as the new HESS data. 



1 Introduction 



The existence of Dark Matter (DM) is by now well established [T], yet its identity is 
a complete mystery; it has no explanation in the framework of the Standard Model of 
particle physics. 

Recently, several exciting data have been reported on cosmic-ray electrons and positrons, 
which may be indirect signatures of decaying/annihilating DM in the present universe. 
The PAMELA collaboration reported that the ratio of positron and electron fluxes in- 
creases at energies of ~ 10-100 GeV [2], which shows an excess above the expectations 
from the secondly positron sources. In addition, the ATIC collaboration reported an ex- 
cess of the total electron + positron flux at energies between 300 and 800 GeV [3]. (cf. 
the PPB-BETS observation [I].) 

More recently, the Fermi-LAT collaboration has just released precise, high-statistics 
data on cosmic-ray electron + positron spectrum from 20 GeV to 1 TeV [5], and the 
HESS collaboration has also reported new data [6]. Although these two new data sets 
show no evidence of the peak reported by ATIC, their spectra still indicate an excess above 
conventional models for the background spectrum [5], and the decaying/annihilating DM 
remains an interesting possibility [7j. 

In the previous studies of DM interpretations of the cosmic-ray electrons/positrons, 
the analyses have been done by (i) assuming a certain source spectrum from annihilat- 
ing/ decaying DM, (ii) solving the propagation and (iii) comparing the predicted electron 
and positron fluxes with the observation. In this letter, we discuss the possibility of 
solving its inverse problem, namely, try to reconstruct the source spectrum from the ob- 
servational data. For that purpose, we adopt an analytical approach, since the inverse 
problem would be very challenging in the numerical approach. 

We show that the inverse problem can indeed be solved analytically under certain 
assumptions and approximations, and provide analytic formulae to reconstruct the source 
spectrum of the electron/positron from the observed flux. As an illustration, we apply 
the obtained formulae to the electron and positron flux above ~ 100 GeV for the just 
released Fermi data [5], together with the new HESS data |6J. It is found that the 
reconstructed spectrum for high energy range is almost independent of the propagation 
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models and whether the DM is decaying or annihilating. We also estimate the errors of 
the reconstructed source spectrum. The obtained result implies that electrons/positrons 
at the source have a broad spectrum ranging from 0(1) TeV to (9(100) GeV, with a broad 
but clear peak around E ~ 400 GeV and another peak at E ~ 1 TeV. 

2 Inverse problem of cosmic-ray electron and 
positron from dark matter 

Let us first summarize the procedure to calculate the e ± fluxes. The electron/positron 
number density per unit kinetic energy f e (E,r,t) evolves as [8] 

§- t fe(E,r,t) = V-[K(E,7)Vf e (E,f,t)} + -^[b(E,7)f e (E,f,t)]+Q(E,f,t), (1) 

where K(E,r) is the diffusion coefficient, b(E,f) is the rate of energy loss, and Q(E,r,t) 
is the source term of the electrons/positrons. The effects of convection, reacceleration, 
and the annihilation in the Galactic disk are neglected. We only consider the electrons 
and positrons from dark matter decay/annihilation, which we assume time-independent 
and spherical; 

Q(E,r-)= q (\f\)^l, (2) 

where dN e /dE is the energy spectrum of the electron and positron from one DM de- 
cay/annihilation, and q(\r\) is given by 

q{\r\) = p(l^1) f° r decaying DM, (3) 

m X T X 



q(\f\) = i^-p(|^) 2 for annihilating DM, (4) 



where mx and tx are the mass and the lifetime of the DM, and (av ) is the average 
DM annihilation cross section. We adopt a stationary two-zone diffusion model with 
cylindrical boundary conditions, with half-height L and a radius R, and spatially constant 
diffusion coefficient K(E) and the energy loss rate b(E) throughout the diffusion zone. 
The diffusion equation is then 

K{E)V 2 f e {E,r) + ^[b(E)f e (E,r)}+ q (\r])^l = 0, (5) 
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with the boundary condition f e (E,r) = for r = \J x 2 + y 2 = R, —L < z < L and 
< r < R, z = ±L. As shown in Appendix A, its solution at the Solar System, 
r = r Q ~ 8.5 kpc, z = is given by 

1 

W)j 



f e (E) = f e (E,r Q ) 



E ^dE' d -^l g (L(E')-L(E)), 



(6) 



where 

dnm 

L(E) 



oo 

£ J ° (i 



n,m=l 

2 



r©\ . /mn\ _ d 
n Rj Sm \~) qnmC ' 



(7) 



Jl(jn)* JO 



dff I" dz J (j n r) sin (|(7r - z)) ? [ J (i?f)2 + j , (8) 



V 71 J 



jl_ m 2 TT 2 
R 2 4L 2 ' 

E dE> K ^ 



b(E') 



(9) 
(10) 



where j n are the successive zeros of Jo- The electron/positron flux is given by & e (E) = 
(c/4n)f e (E). 

Our main purpose is to solve the inverse problem, namely, to reconstruct the source 
spectrum dN e (E)/dE for a given f e (E). As shown in Appendix [Bj this inverse problem 
can indeed be solved, and the solution is given by 



dN e (E) dL{E) 



dE 



dE 



dk 
— < 

2vr 



-ikL(E) J-oo 



f™ dwe~ ikw A( 



w 



L dze lkz g(z) 



where 



A{E) = f e (E)b(E) 



and 

or alternatively, 

dN e (E) 1 
dE 

where the function T(x) is determined from g(x). See Appendix [B] for details. 



A(x = —L(E)) = A(E), 



(ii) 



(12) 



(13) 



(14) 



4 



Several comments are in order. First of all, from Eq. (114j) one can see that the source 
flux at an energy _E src depends only on the observed flux at E ohs > E STC . This may 
be counterintuitive, because the solution to the diffusion equation (Q tells us that the 
observed flux at E — E D ^ S is determined by the source flux at E SIC > E ^, s . It can be 
understood by considering the reconstruction of the source spectrum from the highest 
energy and gradually to the lower energy. 

Secondly, as we will see in the explicit examples, for high energy range the above 
formulae can be approximated by the first term of Eq. ( fl4l) : 

dN e {E) 1 dA(E) 

dE g(0) dE ' [ } 

This is because at higher energy the electrons lose their energy quickly and hence only 
local electrons contribute. Technically, at higher energy g(x) can be approximated as 
g{x) ~ g{0). The approximated formula Eq. (1151) is then directly obtained from Eq. (JS]). 
Note that g(0) is given by the local DM density p & ~ 0.30 GeV/cm 3 ; 

( 1 

Pq for decaying DM, 



9(0) = g(r = r ,z = O) 



m X Tx 



p^——7r for annihilating DM. 

y& 2m 2 x B 



(16) 



3 Examples 



In this section, we show some examples. The energy loss rate is taken as b(E) = E ■ /EqTe, 
with Eq = 1 GeV and r E = I0 16 sec, and the propagation models are parametrized by 
K(E) = K {E/E ) 5 which leads to 

^>-r§(§r < 17 > 

We consider the three benchmark models from Ref. [9], M2, MED, and Ml, which are 
summarized in Table 1. The parameters K and 5 are chosen so that the observed B/C 
ratio is reproduced. 
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Models 


R [kpc] L [kpc] 


5 K Q [kpc 2 /Myr] 


M2 
MED 
Ml 


20 1 
20 4 
20 15 


0.55 0.00595 
0.70 0.0112 
0.46 0.0765 



Table 1: The diffusion model parameters consistent with the observed B/C ratio [9]. 




100 1000 
E (GeV) 

Figure 1: The Fermi data [5] and a fitting curve, together with three different background 
spectra. 



3.1 Reconstruction with the Fermi data 

As the first example, we consider the Fermi data of the the electron + positron flux [5], 
at energy above 100 GeV. As the astrophysical background, we take a simple power 
low oc E~ a , and vary the index as a = 3.2 ±0.1 to represent the effect of the 
background uncertainties^ The normalization is determined by fitting the data for E < 
100 GeV, assuming that the electron + positron flux is dominated by the background 
in this energy ranged For illustration, we fit the observed data by a polynomial as 
E 3 $ e (E) = Y^n=7 c nE n with n max = 5, which is shown in Fig. [1] together with three 



x For instance, for the parametrization of the background spectrum in Ref. |10| . based on the simula- 
tions of Ref. [llj , the high energy electron + positron spectrum is well approximated by a single power 

cfcbg _ ^bg.prim ^bg.sec ^bg.sec ^ ^bg.prim tp-3.25 

2 The PAMELA data [2] shows that the positron fraction is less than O(10) % for E< 100 GeV. 
Assuming that the excess is mainly from the DM decay /annihilation, which generates the same amount 
of electron and positron, at least about 80 % of the total flux is the background in this energy range. 
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Figure 2: The source spectra reconstructed from the Fermi data in Fig. [TJ for different 
background spectrum indices a = 3.1 (top left), 3.2 (top right), and 3.3 (bottom). In each 
case, the results of the analytic formula [ 01 ip or (I14p ] are shown for the decaying DM with 
the three propagation models. The results for the simple approximation formula (1151) are 
also shown, which almost overlap with the Ml model lines in this figure. 



different background spectra (a = 3.1, 3.2, and 3.3). 

In Fig. [21 the source spectra reconstructed by the analytic formula [ (jTTj) or (jl~4"I) ] are 
shown for the decaying DM, with the three background spectra and the three propagation 
models. Here and hereafter, we assume the isothermal DM distribution p(\f\) = p©(r^ + 
r ©)/( r c + l^| 2 ) with r c = 3.5 kpc. We also assume that there is no excess above the 
energy at which the fitting curve crosses the background. Interestingly, in all cases there 
is a clear but broad peak around 400 GeV. There is another increase of the spectrum for 
E ^ 800 GeV, which is due to the last two data points (see Fig. [1]). As can be seen from 
the figure, the reconstructed source spectrum is almost independent of the propagation 
models, except for the energy range E 1z 200 GeV for the M2 model with background 
index a ^3.2. In Fig. [21 we also show the approximated formula (TL5I) . As discussed in 
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Figure 3: The comparison of the source spectra for the decaying and annihilating DM, 
reconstructed from the Fermi data in Fig. [I] with a background spectrum index a = 3.2 
and MED propagation. The results for the simple approximation formula (|15p are also 
shown. 

the previous section, the simple approximation formula (II 5 1) well reproduces the results 
of the full formula. In some cases, the reconstructed source spectrum becomes unphysical 
(negative) in the energy range of E ^ 200 GeV and 700 ^ E ^ 800 GeV. This is because 
of the relatively hard spectrum of the observed flux in this energy range. (The simple 
formula ffl5|) implies that an observed flux much harder than E~ 2 is difficult to obtain 
from the DM decay/annihilation.) 

In order to see the dependence on the DM distribution (or whether DM is decaying or 
annihilating), we show in Fig. [3] the case of annihilating DM compared with the decaying 
DM, for MED propagation with the background spectrum a = 3.2. As expected, the 
reconstructed spectrum is almost independent of whether the DM is decaying or annihi- 
lating. This is also understood from the approximated formula ffT5]) . 

3.2 Reconstruction with the HESS data 

Next, we apply the formula to the newly published HESS data [5J, shown in Fig. HI As 
for the fitting function, we adopt the broken power law described in [6j. We assume the 
same background spectra j^-C 3 - 2 ^- 1 ) as before. The resultant source spectrum is shown 
in Fig. for decaying DM, isothermal DM profile and MED propagation model. One can 
see that the breaking of the power at E ~ 1 TeV results in a peak at that energy in the 
reconstructed source spectrum. 
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Figure 4: The HESS data and a fitting curve [6], together with the three different back- 
ground spectra. 

3.3 Estimation of the uncertainties 

So far, we have neglected the errors in the observed flux. Here, we estimate the effects of 
the uncertainties in the observed flux. As we have seen in the previous two subsections, 
the simple formula (TT5T) is a good approximation in the high energy range E > 200 GeV. 
Therefore, we use it to estimate the errors of the reconstructed spectrum. The derivative 
dA(E)/dE and the corresponding error S[dA(E)/dE] at E — E^ are estimated by fitting 
the sequential three data points at E — (Ei-i, Ei, Ei+i) with a quadratic function. The 
resultant source spectrum is shown in Fig. [6] for the Fermi and HESS data, with the 
three background spectra. Here, we only adopt the statistical errors of the Fermi and 
HESS data [5j [6]. The broad peaks around 400 GeV and around 1 TeV remain even after 
taking into account the errors in the observed flux. Although it is difficult to precisely 
reconstruct the source spectrum because of the limited data as well as the lack of the 
knowledge of the astrophysical background, Fig. [6] implies that the electrons/positrons 
from DM have a broad spectrum, ranging from 0(1) TeV to (9(100) GeV. For instance, a 
direct two-body decay of DM into electron(s), or a three-body decay into electron(s) with 
a smooth matrix element, would have harder spectrum and does not fit the reconstructed 
spectrum well. The obtained result seems to suggest that the source spectrum has a larger 
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Figure 5: The source spectra reconstructed from the HESS data in Fig. HJ for different 
background spectrum indices a = 3.1 (top left), 3.2 (top right), and 3.3 (bottom). In 
each case, the results of the analytic formula [ 01 ip or (1141) ] are shown for the decaying 
DM with the MED propagation models, together with the simple approximation formula 

(USD. 



soft component, like the one from cascade decays. 



4 Discussion 

In this letter, we have discussed the possibility of solving the inverse problem of the cosmic- 
ray electron and positron from the dark matter annihilation/decay. Some simple analytic 
formulae are shown, with which the source spectrum of the electron and positron can be 
reconstructed from the observed flux. It is shown that the reconstructed source spectrum 
at an energy E STC depends only on the observed flux above that energy, E ohs > E STC . 

We also illustrated our approach by applying the obtained formula to the electron + 
positron flux above 100 GeV for the just released Fermi data [5] and the HESS data 
Assuming simple power-law backgrounds, the reconstructed source spectrum indicates 
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Figure 6: The reconstructed spectrum for the Fermi and HESS data, with errors corre- 
sponding to the statistical errors of their data [5j [6], with background spectrum indices 
a = 3.1 (top left), 3.2 (top right), and 3.3 (bottom). See text for details. 



two peaks, a broad one at around 400 GeV and another one at around 1 TeV. 

It is shown that the reconstructed spectrum at high energy is almost independent of 
the propagation models, and whether it is decaying or annihilating. The effect of the 
errors in the observed flux is also discussed. Within the uncertainties, the obtained result 
implies that the electrons/positrons at the source have a broad spectrum ranging from 
(9(1) TeV to (9(100) GeV, with a large soft component, like the one from cascade decays. 

It is difficult to precisely reconstruct the source spectrum at the present stage, because 
of the limited data as well as the lack of the knowledge of the astrophysical background. 
Future measurements, such as the PAMELA data at higher energy range, will allow better 
understandings. There is also a proposed experiment CALET [13], which can measure 
the electron + positron flux up to 10 TeV with a significant statistics, (cf. [2]). We 
expect that our approach will be a useful tool in shedding light on the mystery of DM. 
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A Solution to the diffusion equation 

The diffusion equation (jHJ) can be solved in the following way (cf. Ref. [12J.) Using the 
cylindrical coordinate, f e (E,f) = f e (E,r, z), where r = ' x 2 + y 2 and 

V 2 - — + -— + — (18) 
dr 2 r dr dz 2 ' 



the diffusion equation (pjj) becomes 



K(E) 



d 2 Id d 2 



dr 2 r dr dz 2 
We expand f e as 



f e (E, r,z) + — [b(E)f e (E, r, z)\ + q(r, z)^^ = 0. (19) 



oo _ 

f e (E,r,z) = fnm(E) Jo (jn^) sin {^-(L-z)), (20) 



n,m=l 



fnm \E) 



j2(j ^ J df,f J d * Jotinr)sin (yl 71 "-^)) fe \E,Rr, ^zj , (21) 

where Jo and J\ are the zeroth and first order Bessel functions of the first kind, respec- 
tively, and j n are the successive zeros of Jo- 111 this expansion, the boundary condition 
f e (E,r,z) = for r = R and z = ±L is automatically satisfied. Conversely, any (suffi- 
ciently good) function which satisfies the boundary condition can be expanded as above. 
The orthogonal relations are 

/ dxxJ (j n x)J (j n 'x) = -JxQnjdnnr, (22) 
Jo * 

J da; sin (^-(tt-x)^ sin (^-(tx-x)^ = n5 mm > (m,m'eZ> ). (23) 



Using the differential equation 

i]l ' 1 JoUnx) = 0, (24) 



/ 2 + +3n 

doc oc ditfc 
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one obtains 



-d nm K(E) + ^ + b(E)^ 



„ , . dNJE) . . 

fnm(E) + qnm ^ 1 = 0, (25) 



where d nm and q nm are given by Eqs. and (jHJ). Imposing the boundary condition for 

f nm (E) as 

fnm(E max ) = where £ max = max{£|Q nTO (£) 7^ 0}, (26) 
the solution to Eq. fl25l) is given by 

1 rEmux J AT I rp/\ 

f nm {E) = ^- ) J e dE' ^i-i gnm exp [-d nm (L(E') - L(E))] , (27) 
where £(-E) is given by Eq. ffTUj) . This leads to Eq. ([6]). 

B Solution to the inverse problem 

Here we show two ways of reconstructing the source spectrum dN e (E)/dE. 
By a change of variables, the equation is rewritten as 

A(x) = -J X dy^Mg(x-y), (28) 

where x = —L(E), y = -L(E'), N e (y = —L{E')) = N e (E>), A(E) = f e (E)b(E), and 
A(x = —L(E)) = A(E). We assumed L(E max ) = 0. Here, A(x) and g(x) are known 
functions once we fix the propagation model and have the observational data. Our inverse 
problem is now reduced to the problem of solving for the function dN e (x)/dx, given the 
functions A(x) and g(x). 

The integral equation (|28|) is of the form known as the Volterra's integral equation. It 
is known that the solution of the equation exists and is unique. 

We see that the right-hand side of Eq. fl28l) is a convolution of two functions, so the 
Fourier transform or the Laplace transform seems to be useful tools. In the following, we 
solve the equation in two different ways using the Fourier and Laplace transforms. 
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Fourier transform 

Inserting two step functions in the integral, we get 



A(x) = - dy^Me(y) g{x-y)0(x-y). (29) 

J-oo dy 

Operating dxe~ lkx on each side of the equation and changing the variable as 
z = x — y in the right-hand side, we obtain 



dxe~ tkx A(x) = - / dye~ lky ^' - 9(y) / dze~ lkz g{z)6{z). (30) 

After dividing each side by f^° oo dze~' lkz g(z)9(z) and operate J*^° |^e ifea: , we finally get 
the formula we desire as 

dNe(x) r dk lkx JZdwe-^Ajw) 

or 



v ' J_ 00 2tt J^ o dze-^g(z)9(z) 



V{-L{L)) = — — / — e f oo,.. ih , x • (32) 



d£ v v " d£ J_ 0O 27r f™dze~ ikz g 

Laplace transform 

Eq. fT28l) can be solved in an alternative way. We first differentiate the both side of the 
equation with respect to x and then divide by g(0). Then we obtain 

F(x) = <p(x)- f dyK(x-y)<p(y), (33) 
Jo 



where 



and 



**) = -±>, (34) 



K(x) = -4^-. (36) 
v ' g(0) dx v ; 



The Laplace transform of Eq. (1331) is 

CF(()=£<p(Z)-CK{()C<p(t), (37) 
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where we denote the Laplace transform of a function f(x) by £/(£): 

POO 

£/(0= / dxe~^f(x). (38) 
Jo 

This equation can be solved as 

<p(x) = F(x) + f dyT(x-y)F(y), (39) 
Jo 



where 

rtx) = r £A 



i-ac.™ (40) 

and C~ l denotes the inverse Laplace transform. In the original variable, the solution can 
be written as 

djm = i\d_A(E) + dm i*. dE , d _M Elr _ j 

dE g(0) 
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